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A step by step procedure to derive analytically the exact dynamical evolution equations of the 
probability density functions (PDF) of well known kinetic wealth exchange economic models is 
shown. This technique gives a dynamical insight into the evolution of the PDF, e.g., allowing the 
calculation of its relaxation times. Their equilibrium PDFs can also be calculated by finding its 
stationary solutions. This gives as a result an integro-differential equation, which can be solved 
analytically in some cases and numerically in others. This should provide some guidance into the 
type of probability density functions that can be derived from particular economic agent exchange 
rules, or for that matter, any other kinetic model of gases with particular collision physics. 
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I. INTRODUCTION 

The aim of mechanical statistics is to describe systems 
macroscopically based on the microscopic description of 
interactions of particles within the system. By defining 
the microscopic individual interaction of a collection of 
particles we can describe the system macroscopically by 
means of the Probability Density Function (PDF). A par- 
ticular microscopic interaction will give rise to a definite 
PDF. And vice versa, a given PDF can come only from 
a small set of specific particle interactions. In this re- 
spect, the macroscopic PDF is also providing us with in- 
formation about the microscopic interactions. A classical 
example is the Maxwell-Boltzmann distribution, which 
can be obtained as a solution of the Boltzmann integro- 
differential equation, which he proposed to explain the 
evolution of the PDF for a dilute gas. 

Relatively recently, there has been growing interest in 
reproducing the PDF of money in a real economic system 
by molecular dynamics simulations [T], [2] , [3J , [I] , [5], [5J , 
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[7j. In these statistical models, agents are allowed to ex- 
change money following an exchange rule. PDFs from 
real economies follow Gibb's exponential functions or 
Pareto's laws [5]. Reproducing Gibb's exponential func- 
tions with molecular dynamics simulations has proven 
simple by using a simple money exchange rule that con- 
serves the total quantity of money [5] • Pareto's law dis- 
tributions still remain a challenge, although some results 
do approximate it [9j.[10j. [TIj. p]. 

In this paper, we will show how to derive analytically 
the dynamical evolution equation of the PDF from the 
interaction rules of the agents. This process is very simi- 
lar to the derivation of the Boltzmann integro-differential 
equation from the basic particle collision physics. Let us 
recall that this latter system has a Maxwell-Boltzmann 
PDF in equilibrium. Repetowicz et al. [T3J have de- 
rived, by using mean-field approximation, the first mo- 
ments of some of these economic models based on the 
formal solutions to the nonlinear Boltzmann equations 
derived by Ernst [2]. Lallouache et al. [T5] also calcu- 
late the moments of the random gas-like economic model 
with savings and its steady state solution. During et 
al. derive the moments of some economic models and 
their relaxation times [16]. We will demonstrate a sim- 
ple step by step technique, different from the ones above, 
to derive the integro-differential equations of three eco- 
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Dratnilescu & Yakovenko's model 
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FIG. 1. Dragulescu and Yakovenko's economic model [2] 
which gives an exponential as the steady state PDF. Dots are 
the result of the molecular dynamics numerical simulations, 
and solid line is the exponential function which approximates 
the steady state PDF (Eq. |2]). For these calculations, we take 
< u >= 1. 



money exchange rules in which money is conserved in any 
transaction. 

Assuming that we have N agents trading with each 
other, the mapping that describes this statistical model 
is, 



u'i = e (m + uj), 
u' = (1 - e)(v,i +uj), 
1...N, 



3 

i,3 



(1) 



where e is a random generated number in the interval 
(0, 1), Ui and Uj are the initial money (or energy) and u[ 
and u'j are the final ones of agents i and j, respectively, 
where the pair of agents is also randomly chosen 

for each transaction. These quantities, Ui, will always 
be positive. The steady state distribution, f eq , obtained 
with numerical simulations of this system is shown in 
Fig. [T] (dotted line). It can be easily verified that it is an 
exponential or Gibbs distribution, 



nomic models present in the reviews by Patriarca et al. 
[17] . Chatterjee and Chakrabarti [18 or Yakovenko and 
Rosser [T5]. Some of these systems exhibit exponential 
PDFs as their steady state solutions. Other systems, 
which seem to follow Gamma distributions [SO] at their 
steady state, are in fact only approximated by them and 
do not follow them exactly [13] . We will show that some 
systems, in some particular cases, deviate appreciably 
from these Gamma distributions. We will derive the an- 
alytical formulas to calculate these distributions and we 
will compare these PDFs with Gamma distributions and 
molecular dynamics simulations. This technique should 
shed some light into the relationships between PDFs and 
the microscopic interactions of particles, and in particu- 
lar, into the derivation of a Pareto type distribution from 
molecular dynamics simulations. 

We also present, to our knowledge, for the first time, 
the dynamical equations of the evolution of the PDFs for 
these systems. This will allow us to study these systems 
not only from the point of view of their stationary PDFs, 
but also how they evolve in time. This area of research 
has often been neglected in the past, but it is key if we 
want to model the economy adequately by describing its 
evolution in time and not only in a static way. As an 
example, we will show how to derive the relaxation times 
of these systems. 



II. PURE RANDOM GAS-LIKE ECONOMIC 
MODEL 

We will introduce the technique directly with an ex- 
ample. In this case, we take the mapping introduced by 
Dragulescu and Yakovenko [5] to model the flow and dis- 
tribution of money. This example is one of the simplest 



f eq (u)=/3cxp(-(3u), (2) 

where f(u)du denotes the PDF, or probability of find- 
ing an agent with money (or energy) between u and 
u + du. This PDF is normalized and the mean energy 
per particle then turns out to be 



(«) = jj. (3) 

Let us now derive this function analytically. To do this, 
it is not only important the money exchange mapping 
(Eqs. [lj in the determination of the PDF, but also the 
selection rule to choose which particles do interact at 
each time step. Therefore, we will explicitly describe the 
detailed and complete algorithm, 

1. A pair of different agents are randomly selected in 
the system from a uniform distribution in the in- 
terval [1, A]. This is the pair in Eqs. [I] 

2. A random number e between and 1 is generated 
from a uniform distribution. 

3. Exchange rules of Eqs. [T] arc finally applied. 

4. The application of all of the above "rules" will be 
denoted as a "step". These steps are successively 
and indefinitely repeated. 

We now derive the analytical expression for the PDF 
variation after each step, df(u)/dt, for a given money 
u. Note that with this notation we assume that the time 
unit is one step. As explained in the next paragraphs 
in more detail, it can be straightforwardly seen that this 
variation comes on one side from the probability that 
agent u has to be selected for a particular exchange, i.e., 
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that Ui or Uj take the value u, and from another side 
from the probability that the result after the trading is 



u, i.e. 



that 



or u'j give the value u. We can then write, 



As indicated by Eq. [4j to obtain the final result for a 
generic value u we combine the decreasing term of Eq. [5] 
with the positive contribution of Eq. [8] to get 



df(u) 
dt 



dt 
dt 



dt 
dt 



(4) 



Ui—UA —u. —u. —u 



where all terms are maintained separated in the case 
that the trading rule is not symmetric in the indices (i, j) 
before or after the interaction. 

Let us see the detailed explanation for the present ex- 
ample. In rule [T] one particular agent from the whole 
population N is selected. The probability that this par- 
ticular agent is in the [m, Ui+dui] range is proportional to 
the number of agents in that range, then the PDF for Uj 
will be proportionally depleted by the quantity f(v,i)dUi, 
that is, 



dt 



-fin). 



(5) 



Rule [T] will also deplete the PDF for Uj in a similar 
way, that is, df(uj)/dt ~ —f(uj). 

Rules [2] and [3] imply that, since e is between and 1, 
and according to Eqs. [T] the net result verifies 



N d M = -2f(u) + 2 



dU 



u f( Ul )f(U- Ul ) 

dm , 



(9) 

where N is the normalization factor for df(u)/dt. 
To derive the steady state solution we just set 
df(u)/dt = to finally obtain 



feq(u) 



dU 



, feq{v>l)feq(U - U X ) 

dm . 



(10) 



It is now easy to verify that the normalized exponential 
distribution from Eq. [2] satisfies this last equation. This 
analytical solution is plotted in Fig. [T] where it is verified 
that it is very much in agreement with the numerical 
simulations. 

We shall now derive the relaxation time of this system. 
We will assume that we start with a PDF that is different 
from the equilibrium one (Eq. [2]), but close to it. If we 



10 



substitute this PDF into the right hand side of Eq 
we can verify (not shown here) that we obtain a very 
good approximation to the exponential equilibrium PDF 
(Eq. [2]). The evolution equation close to the equilibrium 
is then simplified to, 



< u\ < Uj + u 



j ' 



(6) 



where u\ is equally distributed in the [0, Ui + Uj] in- 
terval. This implies that the PDF for u\ is increased 
proportionally to the number of agents with U{ and Uj 
and inversely to their total money itj + Uj, that is, to the 
rate f(ui)duif(uj)duj/(ui + Uj). Then, the total varia- 
tion dffu'^/dt will be obtained by integrating amongst 
all the possible values of Ui and Uj that give rise to the 
result u\. In order to find the limits of integration in 
a simple way, we make a change of variables replacing 
(Ui,Uj) by (ui,U), where 



U 



(7) 



Equation|6]is then transformed into < u' t < U, which 
forces the integration limits of U to be between u\ and 
00. If t/ is now fixed, and using its definition from Eq. [7J 
the integration limits of Ui have to be between and U. 
With these new variables, the expression for df^u'^/dt is 
finally obtained 



dt 



du ( u duA u ^ u - u * 

Jo 



u 



(8) 



By symmetry of Eqs. flj a similar result for df(u'j)/dt 
is just obtained by substituting u\ for u'j in Eq. [8] 



N^^^-2f(u) + 2f eq (u). 



(11) 



From here, we can immediately see that the relaxation 
time of the system is 



N 
~2 



(12) 



collisions. 



III. RANDOM GAS-LIKE ECONOMIC MODEL 
WITH SAVING 

We will now explore the gas-like economic model with 
saving introduced by Chakraborti and Chakrabarti pQ. 
The trading rules in this model are given by a random 
mapping that conserves the amount of money in each 
exchange, and in which each agent saves a fixed fraction, 
A, of the money he owns before the transaction. More 
precisely, money is exchanged in the form, 



i,j = 1 



A u 
A u 



e(l 

3-r(l- 

. N. 



- A) (ui + uj), 
e) (1 - A) (ui + Uj 



(13) 
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FIG. 2. Chakraborti and Chakrabarti's economic model [I], 
which gives approximately a Gamma function as the steady 
state PDF, plotted for different values of the saving parame- 
ter, A. Dots are the result of the molecular dynamics numer- 
ical simulations, blue solid line is the Gamma function which 



approximates the steady state PDF (Eq. 14 1 and red dashed 
line is the PDF which is the solution to the integro-differential 
equation |22| For these calculations, we take < u >= 1. 



The steady state distribution obtained by numerical 
simulations of this system is shown in Fig. [2] (dotted 
lines). It is verified |17j that it can be approximated by 
a Gamma distribution, 



feq(u) ~ au n 1 exp {—mi/ (it)) 



(14) 



where (it) is the mean wealth of the multi-agent sys- 
tem, and n is given by 



1 + 2A 
1 - A ' 



(15) 



The factor a is found after the normalization of the 
PDF, 



r(n) 



(16) 



This Gamma function is plotted as a blue solid line in 
Fig. [2] It can be observed that it adjusts relatively well 
to the numerical simulations. 

To derive the analytical expression of the steady state 
PDF we have to look at the precise rules for this mapping: 

1. A pair of different agents of the system are ran- 
domly selected from a uniform distribution in the 
interval [1, AT]. This will be the trading pair 

in Eqs. [T3j 

2. A random number e between and 1 is generated 
from a uniform distribution. 



3. The saving parameter A is maintained constant 
through the whole process. 



4. Exchange rules of Eqs. 13 are finally applied 



5. The application of all of the above "rules" will be 
denoted as a "step". These steps are successively 
and indefinitely repeated. 

Let us derive now in an analytical way the PDF for 
this system. As in the previous model, rule [l] will deplete 
the density function f(u) in by a quantity df{ui)/dt 
given by 



dt 



-/(V 



(17) 



and similarly, df(uj)/dt 

The mapping equations (Eqs. 13 1 lead to the following 
inequalities 



Xui <u[< Xui + (1 - X)U, 



(18) 



where v! i , by using again U = ut + uj , is equally dis- 
tributed in the [Au£,Aitj + (1 — \)U] interval. This im- 
plies that the PDF for u[ is increased proportionally to 
the number of agents with Ui and Uj and inversely to the 
length of the interval, (1 — X)U, where it is distributed, 
that is, to the rate f(ui)duif(uj)duj/((l — X)U). Then, 
the total variation dfiu'^/dt will be obtained by inte- 
grating amongst all the possible values of Ui and Uj able 
to give rise to the result u' { . The limits of integration can 
be easily found from Eq. [l8j 



^ < u'J\, 

Ui > « - (1 - x)u)/x 



(19) 



As in the former section for the Dragulescu and 
Yakovenko model [2], we obtain the result: 



dU 



min[uj /A, U] 

'. — fl — A)C7 1 
>0j 



dti 



g/K) 

dt 

(i-x)u ■ 



(20) 



A similar formula can be derived for the u'j variable. 
Combining these results (Eqs. [l7| and |20|) , the final ex- 
pression for the PDF variation of(uJJot in a generic 
value u is derived: 



N 



df(u) = 

at 

-2f(u) + 



dU 



m'm[u/ X,U] 



f( Ul )f(U- Ul ) 

(l-A)C/ ' (21) 
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where N is again the normalization factor of df(u) /dt. 
To find the steady state solution, we just set 
df (u) /dt = to finally get 



feg(u) 



dU 



i/X,U] 



dui 



feq(ui)feq(U ~ Uj) 

(l-A)t/ 



(22) 

Eq. [22] can be solved iteratively by feeding as first ap- 
proximation the gamma function from Eq. [14] in its right 
hand side and solving the integrals numerically to give a 
second order PDF approximation in the left hand side. 
We need only to iterate once to get a better approxima- 
tion than the first one. Results for different parameters 
of A are shown in Fig. [2] as a dashed red line, together 
with molecular dynamics simulations and the first ap- 
proximation itself (Eq. [14] ). It can be verified how the 
the Gamma distribution deviates substantially from the 
molecular dynamics simulations for the cases with a low 
value of A. On the contrary, the numerical solution of the 
integro-differential equation (Eq. 22 ) matches very well 
the molecular dynamics simulations. A similar expres- 
sion of this steady state solution, derived in a different 
way, was found by Lallouache et al. (Eq. 34 in |15j). 

As before, we can now calculate the relaxation time by 
starting with a PDF close to equilibrium. In this case 



we will start with a Gamma function (Eq. 14 ) , which we 
know is not the distribution at the equilibrium, but it is 
close to it. As we have seen(Fig[2]) , if we place this PDF 
on the right hand side of Eq. [22]we obtain a very good ap- 
proximation of the equilibrium distribution. That leaves 
the evolution equation (Eq. |21[ ) close to equilibrium sim- 
plified to, 
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FIG. 3. Angle-like's economic model [21] which gives approx- 
imately a Gamma function as the steady state PDF plotted 
for different values of uj. Dots are the result of the molec- 
ular dynamics numerical simulations, blue solid line is the 
Gamma function which approximates the steady state PDF 
(Eq. 141, red dashed line is the PDF which is the solution 
to the integro-differential equation |35| For these calculations, 
we take < u >= 1. 



2. We obtain a random number between and 1, e, 
generated from a uniform distribution. 

3. The exchange parameter ui is maintained constant 
through the whole process. 

4. Agent Uj will now give some money to agent u$ 
according to the following mapping, 



N 



df(u) 
dt 



-2f(u)+2f eq (u), 



(23) 



which again leads us to a relaxaion time of r = A^/2 
collisions. 



3 

i,3 



Ui + e uj Uj, 

Uj — € UJ Uj , 

1...N. 



(24) 



IV. ASYMMETRIC RANDOM GAS-LIKE 
ECONOMIC MODEL 

We will now explore a mapping given in |17j as a mod- 
ification of a model introduced by Angle |21j . This map- 
ping has a parameter ui which graduates the amount of 
exchange of money between interacting agents. We will 
only explore the simplest version of all these possible 
Angle-like models, where one agent gives money to an- 
other one regulated by uj and a random number. It is 
an asymmetric model because in this model one of the 
agents is chosen as a winner and the other one as a loser. 
Other more sophisticated versions could as well be ex- 
plored with this method. The precise rules to follow for 
this mapping are, 

1. We select randomly using a uniform distribution 
two different agents 



As before, 



and 



are the amount of 



money before and after the interaction of agent i 
and j respectively. 

5. The application of all of the above "rules" will be 
denoted as "step" . These steps are successively and 
indefinitely repeated. 

The steady state PDF in the numerical simulations of 
this mapping are approximated by a Gamma distribution 
(Eq. ~ 



14 1, where (u) is the mean wealth of the ensemble of 



agents, a is given by Eq. 16 and n verifies the relationship 



3 - 2u) 
2uj 



(25) 



This function, together with a molecular dynamics 
simulation of the model are shown in Fig. [3] Although 



G 



this Gamma function apparently fits very well the nu- 
merical simulation, we will show here that in fact it is 
just an approximation to the PDF and it is only exact 
for the cases where aj = lorw = l/2. 

Let us now obtain the analytical expression for the 
PDF. We will proceed as before. Rule [l] will deplete the 
PDF for those particular values of u 



at J u i uuj 

The final expression (Eq. |4| for the evolution of the 
PDF is 



and 



at 



at 



(26) 



(27) 



From Eqs. [24] we can derive the limits of the variables. 
From the first of these equations ( 24 ) we obtain the in- 
equalities 



Ui < u i < Ui + uj Uj , 



(28) 



from which we see that the variable v! i is spanning an 
interval of length uj uj . This interval is where the prob- 
ability from e spreads over, so we need to divide by this 
factor. We can now split this inequation in the two fol- 
lowing ones 



Ui < u t , 



(29) 



and from the second equation of the mapping ( 24 ) we 
obtain 



Uj > Uj, 



Uj < u'J{l-U)). 



(30) 



The increment of the PDF coming from the first equa- 
tion of the mapping [24] will then be 



9f(K) 

at 



du 



f(Ui)f{Uj) 



f J {u ! - —Ui) I UJ 

and similarly it is obtained 



(31) 



dfju'j, 
dt 



dui 



du 



f(Ui)f(Uj) 



(32) 



where the integration limits of Ui go from to oo be- 
cause there are no constraints on this variable in the sec- 



ond of Eqs. 24 Since /(itj) is normalized we can remove 



it from this expression, 



N 



df(u) 
dt 



2/(«) 
diii 



dit' 



(u— til) / w 

du 2 - 



LOU 2 



LJU 2 



(34) 



where we have again written the dummy variable 
generically as u, u\ and u 2 and N is the normalization 
factor of df(u)/dt. 

The steady state PDF can now be readily derived from 
it by setting df(u)/dt = 0, 



feq(u) = 



2 pit poo 

- dui du 2 

* JO J{u-ux)/uj 

1 [ U/{1 - U) du 2 M^ 

2 Ju uu 2 



feqi.Uijfeqiu-2) 



UJU 2 



(35) 



We can solve this integro-diffcrcntial equation itera- 
tively by substituting f(u) with the known approximate 
solution [14] in the right hand side of Eq. [35] to give a sec- 
ond order approximation in the left hand side. Usually 
one iteration is enough to obtain a more precise solution. 
The solution of this procedure is shown in Fig. |3]as a red 
dashed line. We can see in this figure that this solution 
is indistinguishable from the first approximation in Eq. 
[14] shown in Fig. [3] as a solid blue line, but they are not 
exactly the same, as we will demonstrate now. 

If the Gamma function (Eq. [m] ) is indeed an exact so- 
lution of Eq. [35] the second derivative of the logarithm of 
both expressions should also match. The second deriva- 



tive of the logarithm of the Gamma function (Eq. 14) 



d 2 lnf eq (u) = (n-1) 
du 2 u 2 

If we now multiply by 2u 2 we obtain a constant, 



(36) 



2 u 



2 d 2 ln/ e9 ( M ) 



du 2 



-2(n-l). 



(37) 



Repeating this same calculation in Eq. [35] by introduc- 
ing the Gamma function (Eq. 14) in its right hand side, 
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Angle's model 




are shown in Fig. [4] for different values of uj. As it can 
be seen, the only results that makes this difference equal 
to zero over the whole domain of u are u) = 1 or 1/2. Let 
us prove it. 



If uj = 1 then n = 1/2 (Eq. 25 ) and the first two terms 
of the right hand side of Eq. |39| are zero, leaving Eq. 40 

as 



2(rc- 1) = 0, 



(41) 



which we can readily verify is identically zero. 
If uj = 1/2 then n = 2 and the hypergeometric function 
can be simplified to 



FIG. 4. Plot of 2u 2 



" 2 d dJ ( " ) + 2(n - 1), where the first 
term has been calculated using Eq. |39| as a function of u. 
It is observed that the Gamma function is an exact solution 
only for uj equal 1 or 1/2. For these calculations, we take 
< u >= 1. 



1*1 



■(-l-2<^>) 



r(2n- 1) 



4it 2 / < u > 2 ' 



(42) 



The first term of the right hand side of Eq. 39 is then 



and doing some more elaborate calculations, we find the 
result, 



2u 



2 d 2 In f eq (u) 
du 2 



d Ul (n-l) a[U Ul) . ^ e^^ 1 



e <"> 1 ~ u 



1 



n 



w(l - w)(™- 1 ) V < u > 1 - w 
1 

UJ 



(38) 



The integral can be solved analytically to finally give, 



2u 



2 d 2 ln/ eg (u) = 
du 2 

a(n- l)u n F(n- l)T(n) 
uj n T(2n- 1) 



1F1 \n,2n - 1, 



n uj — 1 
( 

< It > w 



e <"> !-« 
"w(l - cj)(™- 1 ) 
1 



n uj u 



< u > 1 — UJ 



(39) 



where x-Pi is the confluent hypergeometric function. If 
f(u) is really a Gamma function, the subtraction of Eq. 
[39] from the expected result in Eq. 37 should give zero. 
So the question to be answered in this case is: 



2u 



2 d 2 ln/ e? (u) 



du 2 



+ 2(n-l) = 0? 



(40) 



If this result is not zero, it indicates that f eq (u) is not 
exactly a Gamma function. Results of these calculations 



4e _2 < 



1 



< u > 



(43) 



Inserting the values of uj and n, we can readily verify 
that the second term of the right hand side of Eq. [39] is 



4e" 2 <^ 1 + 2 



V < u>J 



(44) 



Adding both together with the last term of Eq. [39] and 
introducing this result in Eq. [40] gives 



4+-+2(n-l)=0, 

id 



(45) 



which we can verify is true. 

To calculate the relaxation times we proceed as before. 
We can simplify the evolution equation (Eq. 34 ) close to 
equilibrium to 



N 



df(u) 
dt 



-2f(u)+2f eq (u). 



(46) 



And again, the relaxation time will be r = N/2 colli- 
sions. 



V. CONCLUSION 

A step by step guide to derive the dynamical evolu- 
tion equations of the PDFs of economic models involving 
money exchange between interacting agents, or for that 
matter, any other type of similar mappings, has been 
explained. The equilibrium distribution can be found 
by exploring its stationary solutions. This leads to an 
integro-differential equation which can be solved analyt- 
ically in some cases and numerically in others. 
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Dragulescu and Yakovenko's mapping [2] can be solved 
analytically giving exactly an exponential distribution 
(see Fig. [lj. Chakraborti and Chakrabarti's model [T] 
can be approximated by a Gamma function, but this re- 
sult is not precise, especially for low values of the sav- 
ings parameter, A. In this case the numerical solution 
of the integro-differential equation provides a better fit 
to the molecular dynamics numerical simulations (Fig. 
[2]). Angle-like's model [5T] steady state PDF can be ap- 
proximated very well by a Gamma function (Fig. [3]) , but 
according to the more exact integro-differential equations 
derived here, it is not an exact solution to the problem 
except for values of the money exchange parameter, lj, 
of 1 or 1/2 (Fig. §). 

This technique should prove useful to get insights into 
the type of PDFs we can expect from a particular map- 
ping, coming this mapping from an economic exchange 
model or from a gas with colliding particles. The ex- 
change rules, in the case of economic models, or the col- 
lision specifications, in the case of gases, will determine 
which kind of PDFs the system will reach in the steady 
state. It is interesting to draw this parallelism, because 
in the economic case, and with the models shown here, 
we have mappings which are not symmetric in time (the 



equations to deplete and increment the PDFs in one time 
step are totally different) which give rise to distributions 
of the Gamma function family. In the case of gases, 
and in particular with the Boltzmann integro-differential 
equation, the collisions are symmetric in time (the equa- 
tions to deplete and increment the PDFs in one time step 
are symmetrical) , giving rise to many interesting proper- 
ties, like a Gaussian function for the steady state PDF 
and the possibility to prove analytically Boltzmann's H 
theorem. 

It should be noted that we are obtaining the time evo- 
lution equations of the PDFs, and not only the stationary 
solutions. This should allow the study of the dynami- 
cal evolution of these systems, something which is key 
to economics, which as we know is not static. This has 
been illustrated in this paper by calculating the relax- 
ation times of these system, but many other applications 
could be devised. 
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